PROGRAM INMSOIL_DRIVER

use ALLOC
use MPI_VARIABLES
use INMSOIL_ARRAYS, only : T, c6, cpool_soil
use INMSOIL_NUMPARAMS, only : ML

implicit none 
include 'soil.inc'

character(len=80) :: fnmap, fnout
integer(4) :: iyear, imonth_old
real(8) :: hour

integer(4) :: coords_syst
integer(4) :: i,j,k
pi = 4.*datan(1.d0)

coords_syst = 2 ! 1 - Cartesian coordinates
                ! 2 - spherical coordinates, the grid is inherited from
                ! ECOCLIMAP

! Grid sizes
nx = 21
ny = 4
ns = 1
nx1 = nx + 1
ny1 = ny + 1
ns1 = ns + 1
ntot=26296
ntime = ntot*30

dx = 9166. ! meters
dy = 9166. ! meters
dt = 3.*3600. !3*3600. ! sec
xlatit = 65. ! latitude of domain senter, degrees
xlongit = 35. ! longitude of domain center, degrees

! Initial time
iyear = 1991
imonth = 6
imonth_old = imonth
iday = 1
hour = 0.

ifsoil = 4
iflake = 0
ifcmodel = 1
radpar = 3
nlakecall = 1
nsoilcall = 1

fcor = 2.*omega*dsin(xlatit*pi/180.)

fnmap = 'map/grid_map'
fnout = 'results'
prt = .false.

size_MPI = 1
rank_MPI = 0
rank_comm3d = 0
dims(1) = 1 ; coords(1) = 0
dims(2) = 1 ; coords(2) = 0
dims(3) = 1 ; coords(3) = 0
horisbound = .true.
lowbound = .true.
upperbound = .true.
xlbound = .true.
xrbound = .true.
ylbound = .true.
yrbound = .true.
xbound  = .true.
ybound  = .true.
call MPA_MYCRDS &        
   & ( &                   
   & nx,    ny,    ns,    &
   & nxsh,  nysh,  nssh,  &
   & nxsh2, nysh2, nssh2, &

   & dims,  coords,       &
 
   & nxi, nxe,            &
   & nyi, nye,            &
   & nsi, nse,            &
   & nx0i, nx1e,          &
   & ny0i, ny1e,          &
   & ns0i, ns1e,          &
   & nx2e, ny2e,          &
   & nxsi, nxse,          &
   & nysi, nyse,          &
   & nssi, nsse,          &
   & nxs2i, nxs2e,        &
   & nys2i, nys2e,        &
   & nss2i, nss2e,        &
   & nx0si, nx1se,        &
   & ny0si, ny1se,        &
   & ns0si, ns1se,        &
   & ns0i_xproc, ns1e_xproc, &
   & ns0i_yproc, ns1e_yproc  &
   & )
nx0i_x = 0
nx1e_x = nx1
ny0i_y = 0
ny1e_y = ny1

call ALLOCVAR_MPI
allocate(ugeos(ns))
allocate(xugeos(ns,5))
allocate(vgeos(ns))
allocate(xvgeos(ns,5))
allocate(sigma0(0:ns+1))
allocate(sigma1(0:ns+1))
ndados = 1
allocate(xrg(0:ndados))
allocate(xra(0:ndados))
allocate(xpl(0:ndados))

! Some surface parameters and initial conditions
call SOLONOII_MPI &
& (xdd2, &
& tsnoi, tsurw, t2noi, &
& wgnoi, w2noi, wrnoi, &
& veg, z0, z0w, z0h, z0hw, &
& xlai, rsmin, alb, zonirr, xlake, watdepth, &
& tsmed, rra, emis, temmin, temmax, evapor, &
& z0hz0, xice, snowini, iclay, isand, iveg, &
& nx, ny, li, lf, fnmap)

nstep = 0
! Time cycle
nmonth_tot=0
kt=0
nrealiz=1
time_cycle : do nstep = 1, ntime

if(kt.eq.ntot) then
kt=1
nrealiz=nrealiz+1
end if
kt=nstep+1208-(nrealiz-1)*ntot

! Updating time parameters
call JULIAN_DATE(iyear,imonth,iday,hour,dt)
  ihour=int(hour)
! Updating surface parameters from ECOCLIMAP when new month begins
if (nstep == 1 .or. imonth /= imonth_old) then
  call READ_ECOCLIMAP &
  & (dx,dy,xlatit,xlongit,imonth,nx,ny,coords_syst,xlake, &
  & watdepth,alb,z0,fcveg,veg,fnmap)
  if(nstep.ne.1)then
  nmonth_tot=nmonth_tot+1
 ! write(*,*) T(2,11,2), potemp(11,2)
  !write the fields to the files
  open (15, file='TSOIL.STD', &
  & access='direct',form='unformatted', recl=ML*(nx1e-nxi+1)*(ny1e-nyi+1), status='UNKNOWN')
  write (15, rec=nmonth_tot)(((T(k,i,j), i=nxi, nx1e), j=nyi, ny1e), k=1,ML)
  close(15)
  open (16, file='C6.STD', &
  & access='direct',form='unformatted', recl=(nx1e-nxi+1)*(ny1e-nyi+1), status='UNKNOWN')
  write(16, rec=nmonth_tot)((real(c6(i,j),4), i=nxi, nx1e), j=nyi, ny1e)
  close(16)
  open (17, file='CPOOL_SOIL.STD', &
  & access='direct',form='unformatted', recl=nsp*(nx1e-nxi+1)*(ny1e-nyi+1), status='UNKNOWN')
  write (17, rec=nmonth_tot)(((real(cpool_soil(k,i,j),4), i=nxi, nx1e), j=nyi, ny1e), k=1,nsp)
  close(17)
  end if

  imonth_old = imonth
endif

!call READ_ERAmet(nstep, nx, ny, nxs2i, nxs2e,        &
 !  & nys2i, nys2e, uwind, vwind,q,potemp,psuf, pl, dt)
uwind(:,:) = 1.d0
vwind(:,:) = -1.d0
q(:,:) = 1.d-3
potemp(:,:) = 290.d0
zlevel(:,:) = 2.d0
psuf(:,:) = 1.d+5

 call READ_ERAmet(kt, nx, ny, nxs2i, nxs2e,        &
   & nys2i, nys2e, uwind, vwind,q,potemp,psuf, pl, dt)
 ! write(*,*) pl(3,3)
   call SURFTO3D &
& (uwind,vwind,q,potemp,zlevel, &
& u,v,qv,pt,pts,phi,phis,pp, &
& phisuf,phig,psuf, &
& rg,rat,pl, &
& tsnoi,t2noi,tsurw, &
& wgnoi,w2noi,wrnoi, &
& xlake,watdepth,tslake, &
& alb,z0w,z0hw, &
& veg,xlai,rsmin, &
& emis,z0,z0h, &
& fcveg, &
& xice, snowini, &
& rra,rn,h,le, &
& gsolo,xdd2, &
& cdm,tsurf,temmin,temmax, &
& evapor,zonirr,uvsuf, &
& deltaz,hsuf, &
& Srad_surf,Lrad_surf, &
& SraddirUV_surf,SraddirPAR_surf,SraddirIR_surf, &
& SraddifUV_surf,SraddifPAR_surf,SraddifIR_surf, &
& fmeth, &
& xugeos,xvgeos,ugeos,vgeos, &
& sigma0, &
& xrg,xra,xpl, &
& ts,t2,wg,w2,wr, &
& dx,dy,xmount,ymount, &
& timeho,dt,dtl, &
& rhoasurf,ptop,xlatit,fcor, &
& iclay,isand,iveg, &
& nx,nx1,ny,ny1,ns,ns1, &
& nstep,ndados,nf, &
& ifsoil,iflake,ifcmodel,radpar, &
& nlakecall, nsoilcall, &
& verbose, &
& li,lf, &
& imonth,iday,ihour,iminu,iseco, &
& fnout,fnmap,prt)

!call READ_ERArad(nstep, nx, ny, nxs2i, nxs2e,        &
!   & nys2i, nys2e, Srad_surf,Lrad_surf,SraddirPAR_surf,SraddirIR_surf,SraddifPAR_surf,SraddifIR_surf,dt)

 Srad_surf(:,:) = 600.
 Lrad_surf(:,:) = 300.
 SraddirUV_surf(:,:) = 0.
 SraddirPAR_surf(:,:) = 100.
 SraddirIR_surf(:,:) = 100.
 SraddifUV_surf(:,:) = 0.
 SraddifPAR_surf(:,:) = 100.
 SraddifIR_surf(:,:) = 100.
 dtl = 2.*dt
   call READ_ERArad(kt, nx, ny, nxs2i, nxs2e,        &
   & nys2i, nys2e,Srad_surf,Lrad_surf,SraddirPAR_surf,SraddirIR_surf,SraddifPAR_surf,SraddifIR_surf,dt)


call SOIL_MPI &
& (u,v,qv,pt,pts,phi,phis,pp, &
& phisuf,phig,psuf, &
& rg,rat,pl, &
& tsnoi,t2noi,tsurw, &
& wgnoi,w2noi,wrnoi, &
& xlake,watdepth,tslake, &
& alb,z0w,z0hw, &
& veg,xlai,rsmin, &
& emis,z0,z0h, &
& fcveg, &
& xice, snowini, &
& rra,rn,h,le, &
& gsolo,xdd2, &
& cdm,tsurf,temmin,temmax, &
& evapor,zonirr,uvsuf, &
& deltaz,hsuf, &
& Srad_surf,Lrad_surf, &
& SraddirUV_surf,SraddirPAR_surf,SraddirIR_surf, &
& SraddifUV_surf,SraddifPAR_surf,SraddifIR_surf, &
& fmeth, &
& xugeos,xvgeos,ugeos,vgeos, &
& sigma0, &
& xrg,xra,xpl, &
& ts,t2,wg,w2,wr, &
& dx,dy,xmount,ymount, &
& timeho,dt,dtl, &
& rhoasurf,ptop,xlatit,fcor, &
& iclay,isand,iveg, &
& nx,nx1,ny,ny1,ns,ns1, &
& nstep,ndados,nf, &
& ifsoil,iflake,ifcmodel,radpar, &
& nlakecall, nsoilcall, &
& verbose, &
& li,lf, &
& imonth,iday,ihour,iminu,iseco, &
& fnout,fnmap,prt)

! write(*,*)Srad_surf(nx/2,ny/2), Lrad_surf(nx/2, ny/2),  tsurf(nx/2,ny/2), h(nx/2,ny/2), le(nx/2,ny/2)

enddo time_cycle

! End of time cycle

write(*,*) 'Reached end of program'

END PROGRAM INMSOIL_DRIVER

